
c   To be called once per fourier mesh level


      subroutine initialise_randoms(iwhichrand,ns1_4,nft)
      implicit none
c-------------- Integer*4 -----------------------------------
      integer*4 ns1_4
c-----------------------------------------------------------
      integer*8 iwhichrand

      integer*8 nft

      if (iwhichrand.eq.1) then
           call setrand(ns1_4)
           return
      endif

      if (iwhichrand.eq.2) then
           call initialise_random_seed(ns1_4,nft)
           return
      endif

      if (iwhichrand.eq.3) then
           call initialise_random_seed_mill(ns1_4,nft)
           return
      endif

      print*,'Impossible value for iwhichrand = ', iwhichrand
      stop

      end


      subroutine derivs(icoord1,icoord2,disp,nft,iwhichrand,lbox,kmin,kmax,
     &   ispheremode, downweight,cic_correct, iquad, local_nx_4, local_x_start_4, 
     &   total_local_size_4)
      implicit none
      include 'mpif.h'
      include 'fftw_f77.inc'
      ! -- rcrain
      include 'nmax_disp.inc'
      ! -------------
      integer*8 icoord1,icoord2

      ! -- rcrain --
      !real*4 disp(*)
      real*4 disp(0:nmaxdisp-1)
      ! -------------

      integer*8 nft,kmin,kmax,ispheremode
      real*8 downweight
      integer*8 cic_correct, iquad
      real*8 lbox 

c-----------   Integer*4 variables -----------------------------------------
      integer*4 local_nx_4, local_x_start_4
      integer*4 local_ny_after_transpose_4
      integer*4 local_y_start_after_transpose_4, total_local_size_4
      integer*4 nx_4,ny_4,nz_4,use_work_4 
c---------------------------------------------------------------------------
      integer*8 iwhichrand
      integer*8 plan1
      real*8 dummy
      integer*8 i
      integer*4 myid,ierr

      call MPI_COMM_RANK(MPI_COMM_WORLD, MYID, IERR)

      use_work_4 = 0
      nx_4 = nft
      ny_4 = nft
      nz_4 = nft

      if (MYID.eq.0) print*,' -in derivs, pre-mpi stage'
      call flush(6)

      call rfftw3d_f77_mpi_create_plan(plan1,MPI_COMM_WORLD,
     &   nx_4,ny_4,nz_4,FFTW_COMPLEX_TO_REAL,FFTW_ESTIMATE+FFTW_IN_PLACE)

      if (MYID.eq.0) print*,' -created plan'
      call flush(6)

      call fftwnd_f77_mpi_local_sizes(plan1, local_nx_4, local_x_start_4,
     & local_ny_after_transpose_4, local_y_start_after_transpose_4,
     & total_local_size_4)

      if (MYID.eq.0) print*,' -set local sizes'
      if (MYID.eq.0) print*,' -total_local_size_4 = ',total_local_size_4
      call flush(6)

      do i=1,total_local_size_4
        disp(i) = 1.e30  ! Crazy values - should all be overwritten
      enddo

      if (MYID.eq.0) then
         print*,'Setting up modes ...'
         call flush(6)
      endif

      call MPI_BARRIER(MPI_COMM_WORLD,ierr)

      call set_up_complex_array(disp,local_x_start_4,local_nx_4,nft,
     & lbox,kmin,kmax,ispheremode,downweight,cic_correct,iquad,iwhichrand,icoord1,icoord2)



      call MPI_BARRIER(MPI_COMM_WORLD,ierr)
      
      if (MYID.eq.0) then     
         print*,'Now doing fft '
         call flush(6)
      endif

      call rfftwnd_f77_mpi(plan1,1,disp,dummy,use_work_4,FFTW_NORMAL_ORDER)

      call rfftwnd_f77_mpi_destroy_plan(plan1)

      call MPI_BARRIER(MPI_COMM_WORLD,ierr)

      if (MYID.eq.0) then  
         print*,'Finished fft'
         call flush(6)
      endif

      return
      end




      subroutine set_up_complex_array(disp_complex,local_x_start_4,local_nx_4,nft,
     & lbox,kmin,kmax,ispheremode,downweight,cic_correct,iquad,iwhichrand,icoord1,icoord2)
      implicit none
      real*8 PI
      parameter (PI=3.1415926535897d0)
      integer*4 local_x_start_4,local_nx_4


      integer*8 nft
      integer*8 kmin,kmax,ispheremode
      real*8 downweight
      integer*8 cic_correct,iquad
      integer*8 iwhichrand, icoord1,icoord2
      real*8 lbox,powspec,kval,kmag
      real*4 ggvar1(0:10240),ggvar2(0:10240)
      complex*8 g_comp
      real*8 pow_fac
      real*8 prefactor
      real*8 fx,fy,fz,arg
      real*8 f1,f2,f3,fh1,fh2,fh3
      real*8 arg2
      real*8 f2xyz   
       

      integer*8 kmag2
      

      integer*8 ii,jj,kk,rcube
      integer*8 lvec(3)
      integer*8 ind1

      complex*8 disp_complex(0:*)

      complex*16 factor1,factor2

      real*8 twopi
      parameter (twopi=2.0d0*3.1415926535897932d0)
      

      !rcrain
      prefactor = 0.0d0

      do ii=local_x_start_4,local_x_start_4+local_nx_4-1
       do jj=0,nft-1
       if (iwhichrand.eq.1) then
          call get_ij_k_rands(nft,ii,jj,ggvar1,ggvar2)
       else
        call get_ij_k_rands23(nft,ii,jj,ggvar1,ggvar2)
       endif

        do kk=0,nft/2
            fx = 1.0d0
            fy = 1.0d0
            fz = 1.0d0
            f2xyz = 1.0d0
            lvec(1) = ii
            lvec(2) = jj
            lvec(3) = kk
            if (lvec(1).gt.nft/2) lvec(1) = lvec(1)-nft
            if (lvec(2).gt.nft/2) lvec(2) = lvec(2)-nft
            if (lvec(3).gt.nft/2) lvec(3) = lvec(3)-nft
            rcube = max(abs(lvec(1)),abs(lvec(2)),abs(lvec(3)))
            kmag2 = lvec(1)*lvec(1) + lvec(2)*lvec(2) +lvec(3)*lvec(3)

            if ((rcube.eq.nft/2).or.(rcube.lt.kmin).or.(rcube.ge.kmax).or.
     &      ((ispheremode.eq.1).and.(kmag2.ge.kmax*kmax)).or.(rcube.eq.0)) then
                  g_comp = cmplx(0.0,0.0)
                  factor1 = dcmplx(0.0d0,0.0d0)
                  factor2 = factor1   
            else
                   if (rcube.eq.kmin) then
                      pow_fac = downweight  
                  else
                      pow_fac = 1.0d0
                  endif

                  if (cic_correct.ne.0) then         ! CIC correction
                    if (lvec(1).ne.0) then
                          arg = PI*dble(lvec(1))/dble(nft)
                          fx  = sin(arg)/arg
                    endif
                    if (lvec(2).ne.0) then
                          arg = PI*dble(lvec(2))/dble(nft)
                          fy  = sin(arg)/arg
                    endif
                    if (lvec(3).ne.0) then
                          arg = PI*dble(lvec(3))/dble(nft)
                          fz  = sin(arg)/arg
                    endif

                          f2xyz = (fx*fy*fz)**2
                  endif

                  if (iquad.ne.0) then !Quadratic interp correction - believed to be exact expression!
 
                    if (lvec(1).ne.0) then
                          arg = PI*dble(lvec(1))/dble(nft)
                          fh1 = dsin(arg)/arg
                          f1  = dsin(2.0d0*arg)/(2.0d0*arg)
                    else
                          fh1 = 1.0d0
                          f1  = 1.0d0
                    endif

                    if (lvec(2).ne.0) then
                          arg = PI*dble(lvec(2))/dble(nft)
                          fh2 = dsin(arg)/arg
                          f2  = dsin(2.0d0*arg)/(2.0d0*arg)
                    else
                          fh2 = 1.0d0
                          f2  = 1.0d0
                    endif 


                    if (lvec(3).ne.0) then
                         arg = PI*dble(lvec(3))/dble(nft)
                          fh3 = dsin(arg)/arg
                          f3  = dsin(2.0d0*arg)/(2.0d0*arg)
                    else
                          fh3 = 1.0d0
                          f3  = 1.0d0
                    endif 

                          f2xyz = (fh1*fh2*fh3)**2 + (fh1*fh1-f1)*fh2*fh2*fh3*fh3 
     &                                             + (fh2*fh2-f2)*fh1*fh1*fh3*fh3
     &                                             + (fh3*fh3-f3)*fh1*fh1*fh2*fh2


                  endif

                  pow_fac = pow_fac/f2xyz    

                  g_comp = cmplx(ggvar1(kk),-ggvar2(kk))
                  kmag = sqrt(dble(kmag2))     
                  kval = (twopi*kmag)/lbox

                  prefactor = pow_fac * sqrt(powspec(kval)/kmag**3/4./twopi)/twopi/twopi/kmag2
                  if (icoord1.eq.0) then
                    factor1 = dcmplx(1.0d0,0.0d0)
                  else
                    factor1 = twopi*dcmplx(0.0d0, lvec(icoord1))
                  endif
                 if (icoord2.eq.0) then

                    factor2 = dcmplx(1.0d0,0.0d0)
                  else
                    factor2 = twopi*dcmplx(0.0d0, lvec(icoord2))
                 endif
            endif

            ind1 = (ii-local_x_start_4)*nft*(nft/2+1) + jj*(nft/2+1) + kk
  
            disp_complex(ind1)   = prefactor* g_comp * factor1 * factor2
 


        enddo
       enddo
      enddo

      return
      end





      subroutine output_psp   ! Write the power spectrum used to file
      implicit none
      real*8 rk,rkval,powspec,powspec_alt

      open (15,file='psp.tex',status='unknown')
       do rk = -3,3,0.01
        rkval = 10**rk
        write (15,*) rk,log10(powspec(rkval)/4./3.1415926/rkval**3),
     &      log10(powspec(rkval))
       enddo
      close(15)

      end



c---------------------------------------------------------------------c
       subroutine fit_spline(name)
       implicit none
       integer*8 npoints
       parameter (npoints=1200)
       real*8 xp(npoints),yp(npoints),yp2(npoints)
       integer*8 np
       common /tran_block/ xp,yp,yp2,np
       integer*8 i
       real*8 yp1,ypn

       character*(*) name

       open (19,file=name,status='old',err=20)

       do i=1,npoints
        read (19,*,end=10) xp(i),yp(i)
       enddo

       stop "Too many data points!"

 10    continue

       np = i-1

       close(19)
       yp1 = 1.e31
       ypn = 1.e31

       call spline(xp,yp,np,yp1,ypn,yp2)


       return
       
 20    continue
       print*,'No transfer function read in '
       print*,'Attempted to read ',name
       stop              ! no spline read in.
       end



       real*8 function powspec(rkval)
       implicit none
       integer*8 npoints
       parameter (npoints=1200)
       real*8 rkval,pow
       real*8 xp(npoints),yp(npoints),yp2(npoints)
       integer*8 np
       real*8 rk 
       common /tran_block/ xp,yp,yp2,np
       include 'powsp.inc'

       rk = log10(rkval)

       if (rk.lt.xp(1)) then
             print*,'rk,xp(1)',rk,xp(1)
             stop 'Out of range!'          ! Large-scales
             return
       endif


       call splint(xp,yp,yp2,np,rk,pow)

       powspec = pow_amp * 10**pow
       return
       end

      subroutine spline(x,y,n,yp1,ypn,y2)
      implicit none
      INTEGER*8 n,NMAX
      REAL*8 yp1,ypn,x(n),y(n),y2(n)
      PARAMETER (NMAX=1200)
      INTEGER*8 i,k
      REAL*8 p,qn,sig,un,u(NMAX)
      if (yp1.gt..99e30) then
        y2(1)=0.
        u(1)=0.
      else
        y2(1)=-0.5
        u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
      endif
      do 11 i=2,n-1
        sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
        p=sig*y2(i-1)+2.
        y2(i)=(sig-1.)/p
        u(i)=(6.*((y(i+1)-y(i))/(x(i+
     *1)-x(i))-(y(i)-y(i-1))/(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*
     *u(i-1))/p
11    continue
      if (ypn.gt..99e30) then
        qn=0.
        un=0.
      else
        qn=0.5
        un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
      endif
      y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)
      do 12 k=n-1,1,-1
        y2(k)=y2(k)*y2(k+1)+u(k)
12    continue
      return
      END


      subroutine splint(xa,ya,y2a,n,x,y)
      implicit none
      INTEGER*8 n
      REAL*8 x,y,xa(n),y2a(n),ya(n)
      INTEGER*8 k,khi,klo
      REAL*8 a,b,h


      klo=1
      khi=n
1     if (khi-klo.gt.1) then
        k=(khi+klo)/2
        if(xa(k).gt.x)then
          khi=k
        else
          klo=k
        endif
      goto 1
      endif
      h=xa(khi)-xa(klo)
      if(h.eq.0.)stop "bad xa input in splint"
      a=(xa(khi)-x)/h
      b=(x-xa(klo))/h
      y=a*ya(klo)+b*ya(khi)+((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**
     *2)/6.

      return
      END



      subroutine normalise_powspec(sigma_8)
      implicit none
      real*8 ll,ul,sum
      include 'powsp.inc'
      external weight_pow
      real*8 weight_pow
      real*8 sigma_8
      real*8 rfilter
      common /radfilt/ rfilter

      pow_amp = 1.0     

      

      ll = -11
      ul =  10

      sum = 0.0


      if (sigma_8.eq.0) then
           pow_amp = 0.0
           return
      endif

      rfilter = 8.0             ! Normalise on 8h^{-1} Mpc scales
      call qtrap(weight_pow,ll,ul,sum)

      pow_amp = sigma_8**2 * pow_amp/sum     ! Normalised variance to sigma_8
c      print*,'pow_amp =  ',pow_amp

c   Check sigma on another scale
c
c      rfilter = 4.0
c
c      call qtrap(weight_pow,ll,ul,sum)
c
c      print*,'sigma_4 = ',sqrt(sum)

      end       

      function top_filter(khimpc)
      implicit none
      real*8 rfilter
      common /radfilt/ rfilter
c   Top hat filter function with  rfilter h^{-1}Mpc sphere

      real*8 top_filter,khimpc, rk

      rk = rfilter*khimpc

      if (rk.eq.0) then
                 top_filter = 1.0
                 return
      endif

      top_filter = 9.*((sin(rk)-rk*cos(rk))/rk/rk/rk)**2
      end



      function weight_pow(logkhimpc)
      implicit none

      real*8 weight_pow,logkhimpc,khimpc,top_filter,powspec

      khimpc = exp(logkhimpc)
      weight_pow = top_filter(khimpc)*powspec(khimpc)

      return
      end



      SUBROUTINE qtrap(func,a,b,s)
      implicit none
      INTEGER*8 JMAX
      REAL*8 a,b,func,s,EPS
      EXTERNAL func
      PARAMETER (EPS=1.e-10, JMAX=25)
CU    USES trapzd
      INTEGER*8 j
      REAL*8 olds
      real*8 dxx,var

      include 'mpif.h'     
      integer*4 myid,ierr

      olds=-1.e30
      do 11 j=1,JMAX
        call trapzd(func,a,b,s,j)
        if (abs(s-olds).lt.EPS*abs(olds)) return
      

      olds=s
11    continue

      call MPI_COMM_RANK(MPI_COMM_WORLD, MYID, IERR)
      print*,'MYID,s,olds ',MYID,s,olds
      print*,'MYID,a,b,j ', MYID,a,b,j
      call flush(6)
      dxx = (b-a)/1000.0d0
      do var=a,b,dxx
       print*,'MYID,var,func',MYID,var,func(var)
      call flush(6)
      enddo
      call system("uname -n")
      call flush(6)


     
      stop 'too many steps in qtrap'
      END



      SUBROUTINE trapzd(func,a,b,s,n)
      implicit none
      INTEGER*8 n
      REAL*8 a,b,s,func
      EXTERNAL func
      INTEGER*8 it,j
      REAL*8 del,sum,tnm,x
      save it

      if (n.eq.1) then
        s=0.5*(b-a)*(func(a)+func(b))
      else
        it=2**(n-2)
        tnm=it
        del=(b-a)/tnm
        x=a+0.5*del
        sum=0.
        do 11 j=1,it
          sum=sum+func(x)
          x=x+del
11      continue
        s=0.5*(s+(b-a)*sum/tnm)
      endif
      return
      END

c-----------------------------------------------------------------
c
c       Alternative powspec routine using linear interp to match L-GenIC
c       Using this routine can improve the level of agreement in velocities
c       with L-GenIC by 3 orders of magnitude for a power spectrum with
c       baryon wiggles. 
c----------------------------------------------------------------

       real*8 function powspec_alt(rkval)
       implicit none
       integer*8 npoints
       parameter (npoints=1200)
       real*8 rkval,pow
       real*8 xp(npoints),yp(npoints),yp2(npoints)
       integer*8 np
       real*8 rk 
       common /tran_block/ xp,yp,yp2,np
       include 'powsp.inc'
       integer*8 j
       real*8 rat


       rk = log10(rkval)

       call locate(xp,np,rk,j)

       if ((j.eq.0).or.(j.eq.np)) then
             print*,'rk,xp(1),xp(np)',rk,xp(1),xp(np)
             stop 'Out of range!'          ! Large-scales
             return
       endif

       rat = (rk-xp(j))/(xp(j+1)-xp(j))
       pow = (1.0d0-rat)*yp(j) + rat*yp(j+1)

       powspec_alt = pow_amp * 10**pow
       return
       end


      SUBROUTINE locate(xx,n,x,j)
      INTEGER*8 j,n
      REAL*8 x,xx(n)
      INTEGER*8 jl,jm,ju
      jl=0
      ju=n+1
10    if(ju-jl.gt.1)then
        jm=(ju+jl)/2
        if((xx(n).gt.xx(1)).eqv.(x.gt.xx(jm)))then
          jl=jm
        else
          ju=jm
        endif
      goto 10
      endif
      j=jl
      return
      END
